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Cosmological Markov Chain Monte Carlo simulation with cmbeasy 

Christian M. Miiller 
Institut fiir Theoretische Physik, Philosophenweg 16, 69120 Heidelberg, Germany 



We introduce a Markov Chain Monte Carlo simulation and data analysis package for the 
cosmological computation package Cmbeasy. We have taken special care in implementing an 
adaptive step algorithm for the Markov Chain Monte Carlo in order to improve convergence. 
Data analysis routines are provided which allow to test models of the Universe against up- 
to-date measurements of the Cosmic Microwave Background, Supernovae la and Large Scale 
Structure. The observational data is provided with the software for convenient usage. The 
package is publicly available as part of the Cmbeasy software at www.cmbeasy.org. 



Introduction 



The wealth of recent precision measurements in cosmolog y I I I I I I I I I I II II ^places 
stringent constraints on any model of the Universe. Typically, such a mode l is given i n te rms 
of a num ber o f cosmological parameters. Numerical tools, such as CmbfastI^, CambI^ and 
Cmbeasy!^, permit to calculate the prediction of a given model for the observational data. 
While these tools are comparatively fast, scanning the parameter space for the most likely 
model and confidence regions can become a matter of time and computing power. The cost of 
evaluating models on a n-dimensional grid in parameter space increases exponentially with the 
number of parameters. In contrast, the Markov Chain Monte Carlo (MCMC) method scales 
roughly linearly with the number of parameters - ■ ■ *. The MCMC method has already been 
used to constrain various rnodels \ A popular tool for setting up MCMC simulations 

is the COSMO-Mc package!^ for the Camb co de, an improved proposal distribution for the 
local Metropolis algorithm has been propose d inl =^. 

We introduce the AnalyzeThis package I2ZI for Cmbeasy. It includes a parallel MCMC 
driver, as well as routines to calculate the likelihood of a model with respect to various data 
sets. We took special care in designing a step-proposal strategy that leads to fast convergence 



^It is part of the cmbeasy v2.0 release. 
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Figure 1: The graphical user interface of Cmbeasy. It can be used to marginalize, visualize and print the one and 
two dimensional likelihoods from the MCMC chains. Shown is the marginalized likelihood in the Q,mh^ — Q,bh^ 

plane of a ACDM model. 



and mixing of the chains. The automated step optimization considerably improves performance 
and is rather convenient. The raw data files can be processed from within a graphical user 
interface (gui). Using the gui, one can marginalize, visualize and print one and two dimensional 
likelihood surfaces (see figure P). 



2 Markov Chain Monte Carlo simulation 

In the following, we will assume the reader is fami liar with the basic ideas of Markov Chain Monte 
Carlo simulationl^^ and the Metropolis algorithnlp^ . We will now describe our implementation 
as contained in the released version!^. 

The global Metropolis algorithm chooses new steps for a Markov Chain via a symmetric 
proposal distribution q{6i-i,9i) = q{9i,6i-i)^ where 9i is a parameter vector that specifies a 
given model. In our implementation, we assume fiat priors P{0) on the parameters, and we 
assign likelihood zero to any parameter set that has at least one point outside the prior. 

A very important aspect of MCMC is to test when the chains are sampling from the under- 
lying distribution. Since at the beginning, the chain migrates from its random starting point 
to regions of higher likelihood, there is a "burn-in" associated with each chain that must be 
eleminated prior to parameter estimation. In principle, it may be difficult to tell from a single 
chain if it has converged. In MCMC, one therefore uses several chains with random starting 
points and monitors mixing and convergence. Our implementation employs the convergence 
test of Gelman and RubiiP^. The key ingredient for this test is a parameter R which can be 



computed from previous chain points. This parameter is a comparison of the variance within the 
chains compared to the variance between different chains. A value of i? < 1.2 for each parameter 
indicates the chains have converged and all previous points should be removed. If one uses the 
gui for chain analysis, the burn-in is automatically removed. 

Since there is no generally accepted procedure to determine when one has generated enough 
chain points for reliable estimation, the algorithm just runs indefinetly in our implementation. 
However, any "breaking-criterion'^^^Jj j^ay be implemented easily. The chains may be monitored 
with external programs during the run. 

The number of steps needed for good convergence and mixing depends strongly on the step 
proposal distribution. If the proposed steps are too large, the algorithm will frequently reject 
steps, giving slow convergence of the chain. If, on the other hand, the proposed steps are too 
small, it will take a long time for the chain to explore the likelihood surface, resulting in slow 
mixing. In the ideal case the proposal distribution should be as close to the posterior distribution 
as possible - which unfortunately is not known a priori. While a simple Gaussian proposal 
distribution with step sizes a^ is sufficient, it is not optimal in terms of computing costs if 
cosmological parameters are degenerate. Instead of using a naive Gaussian proposal distribution, 
we sample from a multivariate Gaussian distribution with covariance matrix estimated from the 
previous points in the chains. By taking into account the covariances among the parameters, we 
effectively approximate the likelihood contour in extent and orientation - the Gaussian samples 
are taken along the principal axis of the likelihood contour. 

The convergence can be further improved by scaling the covariance matrix with a variable 
factor a. Using q, we c an cope better in situations where the projected likelihood takes on 
banana shapes such as in'=^. It also improves the convergence during the early stages when the 
low number of points available limits the estimate of the covariance matrix. We dynamically 
increase a if a chain takes steps too often, while we decrease a if the acceptance rate is too low. 
By this procedure the convergence is speeded up by a factor of about four compared to naive 
Gaussian sampling. 

One can show that modifying the proposal distribution based on previous chain data during 
the run may lead to a wrong stationary distribution. Therefore, we only apply the dynamical 
strategy of finding an optimal step proposal during the early stages of the simulation. When 
the convergence is better than R = 1.2 and the chain has calculated a certain number of points, 
we freeze in the step proposal distribution. 

3 The Software 

The package is part of Cmbeasy and consists of two main components. The first one is a MCMC 
driver using LAM/MPI for parallel execution of each chain. The second one is the AnalyzeThis 
class which is designed to evaluate the likelihood of a given mod el w ith respect to various data 
sets. The se sets in clude t he lat est data of WMAP TT and TeI^, ACBArSI, CBI^I, VSA^l 
, 2dFGRsEE2Iini^ SDSstHKS^ the SNe la compilations of Riess et al. 1^, Tonry et al.El and 
Knop et. al.l^as well as the IfA Deep Survey SNe la data^^. Data files for all experiments are 
included for convenience. New data sets are added continuously to the code. 

The example MCMC driver consists of two routines: master () and slave(). Using LAM/MPI 
for parallel computing, one master and up to ten slaves may be started. The master() will 
determine the initial random starting position for each chain. In a never ending loop, it then 
sends the parameters to the slave()'s and collects the results when the computation is finished. 
Whenever a step has been successful, it stores the parameters and likelihoods of the last step 
together with the number of times the chain remained at the same point in parameter space in 
a new line of one file per chain. The master () monitors convergence and mixing and determines 
the next step for the slave(). Before freeze in, the covariance is estimated and the step proposal 



accordingly modified. After freeze in, the proposal distribution remains unchanged. 

The gui may be used to process the raw output files of the MCMC chains. It enables quick 
and convenient analysis of MCMC simulations. To get started immediately, we include raw data 
from an example MCMC run in the resources directory of cmbeasy. Two and one dimensional 
marginalized likelihoods may then be plotted and printed from within the gui (see figure^. 

4 Conclusions 

We have introduced the AnalyzeThis package, which can be used to constrain cosmological 
models using observational data sets. The AnalyzeThis class provides many functions to com- 
pute the likelihood of a model with respect to measurements of the CMB, SNe la and Large 
Scale Structure. 

In order to constrain models of the Universe with a substantial number of parameters, we 
include a Markov Chain Monte Carlo driver. As the MCMC step strategy determines the 
convergence speed of the chains, we implemented a multivariate Gaussian sampler with an 
additional dynamical scaling. 
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